On the stability of some exact solutions to the 
generalized convection-reaction-diffusion equation 

V. Vladimirov and Cz. Mg,czka 

Faculty of Applied Mathematics, 
AGH University of Science and Technology, 
Mickiewicz Avenue 30, 30-059 Krakow, PL 
O ' Email address: vsevolod.vladimirov@gmai.com 

(N 
u . 

CIh, Abstract Stability of a set of travelling wave solutions to the hyperbolic generalization of the 

' convection-reaction-diffusion equation is studied by means of the qualitative methods and numerical 

00 ' simulation. 

'• 1 Introduction 

d 

In recent decades significant attention was paid to the study of the family of convection- 
reaction-diffusion equations 

ut= [k{u)u^]^ + a{u)u^ + f{u). (1) 

Cf^ ! Equations belonging to this family describe a number of natural phenomena, such as 

^ I transport in porous media, or the motion of a thin sheet of viscous liquid over the 

^ inclined plate (see and the literature therein). This class also contains a nonlinear 

Q ■ generalization of the Focker-Plank equation [2] and a number of models encountered in 

O . the biological sciences [3]. Another source of inspiration for studying the convection- 

reaction-diffusion equations results from the fact that the equation ([T]) represents one of 
the simplest nonlinear models describing phenomena of patterns formation and evolu- 
^ ' tion ^ [5] . It is, perhaps, the combination of relative simplicity and richness of physical 

■ contents, that made the family ([1]) the objective of numerous studies within the sym- 

metry approach, purposed at constructing nontrivial exact solutions and finding out 
the conservation laws [5]-pi]]. 

In this paper, we consider the following evolutionary equation (referred to as GBE): 

autt + ut + fiuux- K,u^x = fiu). (2) 

Here /x, k are positive constants, a is nonnegative, f{u) is a smooth (polynomial) 
function, which will be specified later on. Equation ([2]) is a hyperbolic generalization 
of the convection-reaction-diffusion equation. Let us note, that the term a Uu appears 
when the memory effects are taken into account Equation ([21), as well as its 

numerous modifications, were intensely studied in recent years within the generalized 
symmetry approach [10], |ll]-[l9]. Owing to these studies, the analytical description 
of a large variety of traveling wave (TW) solutions is actually available, including 
interacting traveling fronts, soliton-like solutions, periodic waves, compactons, shock 
fronts and many other. Undoubtedly, knowledge of exact solutions to a non-linear PDE 
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is a great advantage. At the same time, individual exact solution is interesting and 
important from the point of view of applications, if it is typical in some sense to the 
equation under consideration. In most noteworthy cases, self-similar exact solutions 
serve as the intermediate (or the true) asymptotics ^20j-[22j, manifesting attracting 
features. 

The first stage towards the estimation of validity of the exact solution is a study 
of its stability, and this is the main topic of the present work. We formulate the 
conditions which guarantee the stability of some class of TW solutions to the equation 
([2]), obtained in [16j. The structure of the study is following. In section 2 we present 
a family of exact TW solutions, satisfying under certain conditions equation ([2]) and 
formulate the conditions that guarantee the stability of some exact solutions in explicit 
form. In section 3 we construct the numerical scheme based on the Godunov method 
[23]- [25] and bring the results of numerical simulations, backing the qualitative study 
and partly completing it. In the last section we briefly summarize the results obtained 
and outline the ways of further investigations. 



2 Stability analysis of the exact solutions to the 
equation (M) 



2.1 Statement of the problem 

Let us reformulate the results obtained in pjj for the equation ([2]), assuming that 

f{u) = u {u — mi) {u — 7712) {u — 777,3) ; 
where 777^, k = 1, 2, 3, are constant parameters. We are looking for the TW solutions 

u{t, x) = U{z) = U{x-Vt), (3) 

where ^ is a constant velocity of the wave pack. After the formal substitution of the 
travelling wave ansatz ([3]) into the equation ([2]), one obtains a nonlinear second order 
ODE 

{aV^ - k) U + U {fiU ~V) = u {U - mi) {U - 7722) {U - m^) , 

which is, generally speaking, non-integrable. In order to obtain the exact TW solutions, 
we employ a Hirota-like ansatz u{t, x) 
to the following third-order ODE: 



which, being substituted to ([2]), leads 



^2 



A -V"^" - u '^rrii mj"^' + i/ rrii 7772 rrig vE^ 



1=1 



+ 



+ 



+ {2A-fi-u) {^f = 



(4) 
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where IS. = aV"^ — k. For physical reasons [27], we assume that A > 0. 

On first sight, the equation (j4]) is even more comphcated than that obtained by the 
convenient TW ansatz (IHl). But as it easily seen, (jlj) reduces to the linear equation 

A — V — U [{mi 7712 + fTT-l ^3 + ^2 ^3) ^' — ^1 ^2 ^3 ^] = 0, (5) 

provided that 

/i = 3A, (6) 

u = -A, (7) 

V = [mi + 1712 + ^3) A. (8) 

As it was shown in [16j, the roots of the characteristic equation, corresponding to the 
linear equation (Q, coincide with the parameters mi, m2, when the restrictions 
([n])-® take place. This enables us to formulate the following result: 

Theorem 1. Let the equalities ({6])- ([8]) take place. Then, depending on the values 
of the parameters {mfcj^^i, the equation ([2]) has the following exact solutions: 

1. 

_ '"^1 ^1 ["^1 ^] + "^2 C2 exp [m2 z] + ma C3 exp [ms 
Ci exp [mi z] + C2 exp [m2 z] + C3 exp [ms z] 
if mi 7^ m2 7^ m3 7^ mi; 

2. 

_ mi Ci exp [mi z] + exp [m2 ^] [C2 m2 + C3 + ^2 C3 ^] ^^^^ 
Ci exp [mi 2;] + exp [m2 z] [C2 + z C3] 

if mi 7^ m2 = rris] 

3. 

""■"'='"+ C3 + .(C. + ^) 

if mi = m2 = m3 = m; 

4. 

m3 C3 exp [m3 z] +2 exp [a ^] [a cos {/3 z) — /3 sin /3 z] 
C3 exp [m3 2;] + 2 exp [a z] cos (/3 2;) 

if ms is real, while rn2 = mi = a + i /3, a, (3 E R. 



Remark. Usi7ig aTid the iTiequality A > 0, one easily gets the following 

expression for the wave pack velocity: 



1 + Jl + 4aK (^^_^m,) 



2 a ^ mi 
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In order to study the stability of TV solutions, depending in fact on a single variable 
z = t — l^x, it is instructive to pass to new independent variables 

i = t, z = X — V t, 

in which the invariant solutions (P])-f ll2p become stationary. In the new variables the 
equation (|2]) reads as follows: 



a 



dt d z 



u + 



dt d z 



du d"^ u 



u + fiu -— — K 



d z d z"^ 



f{u) (14) 



(for simplicity, we omit the bars over the independent variables henceforth). 

On studying the stability of stationary solutions, we proceed in the standard way, 
presenting the perturbed solution in the form 



u{t, x) = U{z) + e exp [—A t] g{z), 



(15) 



where U{z) denotes one of the TW solutions described by the theorem 1. Up to O(e^), 
the function g{z) satisfies the equation 



A g"{z) + (/i U{z) + 2aVX-V) g\z) + 



a 



TTii rrij 



2 u U{z) J2^i-^^U\z) + n U'{z) 



i=l 



For technical reasons, it is instructive to get rid of the terms proportional to g'{z), and 
this can be done by the substitution 

g{z) = h{z) exp[(p{z)]. (16) 

One easily verifies by the direct inspection, that the following statement holds true: 
Lemma 1. If 

f'{x) ■■ 



lxU{z) + 2aV X - V 



2 A 



then the function h{z) satisfies the equation 



L[h{z),\] = 

Aan h{z) - 4 B{z) h{z) A + h{z) K{z) - 4 A^ h"{z) = 0, 



(17) 



where 



B{z) = K-3aA^U{z)J2 



rrii 



i=l 



K{z) = - 2 ^ mimj + 2U{z)^mi-'iU'^{z) -QU'{z). (19) 



i=l 



i=l 



So, using the ansatz ( !T5|) . followed by the substitution (fT6|) . we get the generalized 
eigenvalue problem (fT7|) . Evidently, the stability of the self-similar solution U{z) can 
be achieved, if all possible values of the parameter A are positive. 

In what follows, we'll restrict our consideration to a family of perturbations, van- 
ishing beyond some compact set < —L, L >. With such a restriction, we get the 
eigenvalue problem 

L[h{z),X] = 0, h{-L) = h{L) = 0. (20) 

Let us note, that in the parabolic case, i.e., when a = 0, fl20l) reduces to the standard 
Sturm-Liouville Boundary Value Problem. Although the eigenvalue problem we deal 
with differs from the classical one, the main conclusions concerning the properties of 
the eigenvectors and eigenfunctions remain the same under quite general assumptions 
[26] . satisfied with certain by the functions B{z) and K{z). 

In order to obtain the restrictions on the signs of the eigenvalues A, we multiply the 
equation (ITTl) by h{z) and then integrate the resulting equation over z from —L to L. 
As a result, we obtain the quadratic equation with respect to A: 



where 



A2-6A + r = 0, (21) 



j':^B{z)h{zfdz A2 j':^K{z)h{zfdz + A^^\\h'\ 



aK\\h\r 4a K, 



J^^h{zy d z, Wh'W^ = J^^h'{zy d z. From the above formulae, we get the 
following relations concerning the roots of the quadratic equations: 

B(z)h(z)^dz 

Ai + A2 = , (22) 

a K 



2 J_^K{z)h{zYdz + 4:\\h'\ 



AiA2 = A2 ^-" y , (23) 

This immediately leads us to the statement: 

Proposition 1. In order that the eigenvalues A^, A; = 1, 2 be positive, it is 
sufficient that the functions B{z) and K{z), restricted to the segment < —L, L >, 
satisfy the following inequalities: 

B{z) > 0, K{z) > 0. (24) 



Below we pose the conditions that guarantee the fulfillment of the inequalities fl24j) for 
some exact invariant solutions to the equation ([2]). 
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2.2 Stability analysis of the solution ([9]) 

We restrict our consideration to the real constants {mi}^^-^^. Without the loss of gener- 
ality, we can assume that they are ordered as follows: < mi < m2 < m^. We assume 
in addition that the constant Ci is nonzero, and the solution ([9]) can be rewritten in 
the form 

U (z) = , "^{z) = exp [mi z] + C2 exp [m2 z] + C3 exp [1713 z] (25) 

The solution occurs to possess the following property: 

Lemma 2. Function (125!) is monotonic for any positive C2 and C3 and satisfies 
the inequalities 

mi < U{z) < 1713. (26) 



Proof. Since the derivative of U (z) is expressed by the formula 

^"(z)^(z) - ^'(z)2 



U'(z) 



we concentrate upon the estimation of the sign of the numerator. After some algebraic 
manipulation, performed with the help of Mathematica package (and which can be 
easily verified manually), we get the inequality 

\E'"(z)\&(2;) - "^'{zY = C2 exp [mg z] {exp [mi z] (mi - mg)^ + 
+C3 exp [m3 z] (m2 - m3)^} + C3 exp [(mi + ms) z] (mi - m3)^ > 0. 

The validity of the inequalities ( l26l) appear from the calculation of limits: 

lim U{z) = 1713, lim U{z) = mi. 

2—^+00 2— >~oo 

The above lemma can be used for the estimation of the signs of inequalities flM|) . 
We begin with the first one. The validity of the inequality B{z) > appears from the 
inequality 

n > a fiV max^g/jf/ (z) = a fiV > a fiV U{z). 

Taking into account the conditions ([6])-([8]), and expressing V by means of the formula 
( IT3l) we can rewrite the inequality k, > a fiV in the form 

r 1^ 

1 + ^Jl + Aa K [mi + m2 + m3)^ 
K > 3 a m,3- — 



4 {mi + m2 + m3)'^ 
or, what is the same, 

Aan [mi + m2 + m3)^ > 
> 3 m3 |2 + 4 a K (mi + m2 + m3)^ + 2 ^Jl-\- Aan [mi + m2 + m3)^ | . 



This, in turn, is equivalent to 

4 a K (mi + m2 + ms)^ [(mi — ms) + (m2 — ms)] > 

> 6m3 |l + a/1 + 4aK(mi + m2 + m3)2| . (27) 

It is evident, that under the above assumptions the inequahty (^7^ cannot be fulfilled, 
so, following this way we cannot gain any useful information. It occurs to be possible, 
if we restrict the set of exact solutions described by the formula (125|) by putting C3 = 0: 

^ _ mi exp [mi z] + C2 exp [ma z] 
exp [mi z] + C2 exp [m2 z] 

In analogy with the lemma 2, one can check the validity of the following statement: 

Lemma 3. if C2 > 0, then the function ( 128|) is monotonic and satisfies the inequal- 
ities 

mi < U{z) < 1712- 

So, now there is the condition k, — a fiV m2 > 0, which guarantees the validity of the 
inequality i?(z) > for the exact solution f l2S]) . Using (EIl-dHD, and their consequence 
(IT^. we get the inequality 

4 a K (mi + m2 + m3)^ [(mi — m2) + (m3 — m2)] > 

> 6 ma 1 1 + a/1 + 4 a K (mi + m2 + rris)'^^ . (29) 

Obviously, the inequality does not hold for arbitrary values of the parameters. Yet 
if all the parameters, but m3 are fixed, then the LHS behaves as mg while the RHS as 
mg. Hence there exists the critical value such that for any ms > mg^ the inequality 
( |29|) does take place. 

Now let us consider the condition 

3 3 
K{z) = ^ m^ - 2 ^ mi mj + 2 U{z) ^ m, - 3 t/^(z) - 6 U'{z) > 0. 

1=1 ij^j i=l 

Below we shall use the following elementary statement: 

Lemma 4. The derivative of the function fl28l) satisfies the inequality 

< U'iz) < (-H^l^^. 

^ ' - 4 

Using the above restrictions and the statements of the lemma 3, we get the estimation 

3 

K{z) > "^ml - 2^^771^771^ + (30) 

fe=l iy^j 

2 (^2 — mi)^ 
rrik — 3 m2 — 6 = Ki. 

k=l 



Fixing all the parameters but ma, we can treat Ki as the quadratic function: 

Ki = ml — 2 7712 f^l + C (mi, 7712) . 

So there exists a number tti^^ such that for 1713 > m*^ , Ki is positive. From this 
immediately follows the main result of this section: 

Theorem 2. If ms > maxjmg^, m*^}, then, under the restrictions stated above, 
the TW solution (EHD is stable. 



3 Numerical study of the invariant TW solutions 
to the equation ([2]) 

3.1 Construction of the numerical scheme. 



We base our numerical calculations on the Godunov method [231 El]. K is not difficult 
to extend the construction of a numerical scheme upon somewhat more general equation 



auu + ut + fiuu^- K = /(m). 



(31) 



containing nonlinear diffusive term. Introducing the new variable = Ut—y/K'yu" Ux, 7 
a~^, we can rewrite the equation (IHT]) in the form of the first order system: 



9_ 



u 



— a/k7 

■ nn/ 2-1 







d_ 

dt 



u 



H, (32) 



where H = {\E', 7 [f{u) — \E']}*'^ , and (■)*'' stands for the operation of transposition. 

Let us consider the calculating cell abed (see Fig. 1) lying between m — th and 
(m + 1) — th temporal layers of the uniform rectangular mesh. It is easy to see that 
the system ( 13 2 p can be presented in the following vector form: 



OF dG 

+ 



dt 



d X 



with F 



{u, M/) 
G = 



/7K- 



u 



n/2+1 



1 
2 



n/2 + 1 

From fl33l) arises the equality of integrals 

'dF dC 
dt dx 




n/2 + 1 



dx dt 



(33) 



+ ^ a/7 k 



tr 




Hdx dt. 



where fl is identified with the rectangle abed. Due to the Gauss-Ostrogradsky theorem, 
integral in the LHS can be presented in the form 




OF dG\ , , 
-rrr + -rr- ' dx dt 

Ot ox 



Gdt-Fdx. 



(34) 



m+l 
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Figure 1: Scheme of the calculating cell 




Let us denote the distance between the i — th and {i + 1) — th nodes of the OX axis 
by Ax while the distance between the two adjacent temporal layers by At. Then, up 
to O (|Ax|^ , |At|^), we get from the equation fl33l) the following difference scheme: 

where i , G"^ i are the values of the vector-function G on the segments b d and 

«+2 *-2 

a c, correspondingly. In the Godunov method these values are defined by solving the 
Riemann problem. Below we describe the procedure of their calculation. 

In accordance with common practice, instead of dealing with the initial system fl5^ . 
we look for the solution of the Riemann problem '^i) at x < and (m2, ^'2) at 
X > to corresponding homogeneous system 



d [ u\ f -Ci \ d 



u 



where 



0, (36) 



Ci = a/7k<. 



, ^^, r- n/2-l , 1/2 3/2 n/2 

C2 = 7 Mo + ^0 V7 2 "0 +k/7^Mo . 
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Using the linearized system (136|) . it is easy to calculate the Riemann invariants 



corresponding to the characteristic velocities C± = ±Ci. Characteristics x = ±Cit 
divide the half-plane t >0 into three sectors (see Fig. |2]) and the problem is to find the 
values of the parameters in sector II, basing at the values (ui, ^'i) and (m2, "^2), which 
are assumed to be defined. The scheme of calculating the values un, "^ii is based on 
the property of the Riemann invariants to retain their values along the corresponding 
characteristics. From this we get the system of algebraic equations 

U2 = Un- 

So the values of the parameters u, \l/ in the sector — Ci t < x < Cit are given by the 
formulae: 



Un = U2, 
■q/jj = ^^J, + C2 ^TTT^. 

Thus, the difference scheme for (1311) takes the following form: 

At 



u 



m+l 



m+l 



Ax 
At 

Ax 



{G2 



{G2 



At {H2 



where 



/-IK 



u: 



n/2+1 



m 

2 



7 
2 



2 2^7^ / 



n/2 + 1 

n/2+1" 



2,3, 



7/? M 



A^- 1, 



n/2 



(37) 



2,3, ...A^ — 1, are calculated by means of the formula (|3 



in which (mi, '^i) and (m2, ^['2) are substituted, correspondingly, by (m^i, and 



\E'™'), while the constants C^, k = 1,2 take the form 



1 = ^7>< (^^1)" . 
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Figure 3: Numerical solution of the system ([2]) in case when the invariant kink- like 
solution (128|) with a = k, = 1, mi = 0.5, m2 = 1.5, ms = 5, Ci = C2 = 1, C3 = 
is taken as the Cauchy data. Successive graphs of TW solution, moving from left to 
right, correspond to = 4 (i — 1), i = I, ...7 

3.2 Results of numerical simulation 

Below we present the results of numerical solution of the Cauchy problem for system 
(|2]). In all numerical experiments the parameters a and k, were taken to be equal to 
one, while the remaining parameters varied from one case to another. In the first series 
of the numerical experiments we put mi = 0.5, m2 = 1.5, and, since for this choice 
777,3 ~ took 7773 = 5 in order to satisfy the requirements of the theorem 2. As 

the Cauchy data we used the invariant solution described by the formula f l28|) . and 
corresponding to to = 0. Results of the numerical simulation are shown in Fig. |3l It 
is seen that the kink-like solution evolves for a long time in a stable self-similar mode. 
Figure H] shows the graphs of the functions B{z) and K{z) for the above values of 
the parameters. It can be seen on this graphs that both of the functions are strictly 
positive in a vicinity of the front of the kink-like solution (!28|l . 

Next we performed the numerical experiments in which the full solution (Q was 
taken as the Cauchy data. We put in this case C3 = 3 and the rest of the parameters 
remained the same as in the previous case. The results of numerical simulation show 
that the self-similar solution is unstable. Fig. El The source of the instability is seen 
in Fig. ini showing the graphs of the functions B{z) and K{z). Both of these functions 
are negative in some vicinity of the origin. Besides, the function K{z) has a local 
minimum in this vicinity, in which it attains a sufficiently large negative value. Analysis 
of the formula ( |T9|) shows that the presence of such maximum can be attributed to the 
abruptness of the slope of the kink-like solution. Note that the graphs of the functions 
shown on Fig. |5] correspond to sufficiently large times tj > 12, for which the effects of 
instability become evident. Therefore they do not coincide with initial profile described 
by the formula ([9]), which is quite sharp. 

Let us briefly describe the results of some other numerical experiments. Fig. [7] shows 
the results of the numerical evolution of a step-like initial perturbation described by 
the formula (Q. Fig. |8] shows the results of the numerical evolution of the initial 
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B(black);K (dashed) 



Figure 4: Graphs of the functions B{z) and K{z), obtained for a = k = 1, mi 
0.5, m2 = 1.5, ms = 5, Ci = C2 = 1, and C3 = 




Figure 5: Numerical solution of the system ([2]) in case when the invariant kink-like 
solution ([9]) with a = k = 1, mi = 0.5, m2 = 1.5, m^ = 5, Ci = C2 = 1, C3 = 3 is 
taken as the Cauchy data. Successive graphs correspond to ti = 12+4 (i — l), = 1, •••4 



B(black);K (dashed) 
15 r 




Figure 6: Graphs of the functions B{z) and K{z), obtained for a = k = 1, mi = 
0.5, 777.2 = 1-5, 7ri3 = 5, Ci = C2 = 1, and C3 = 3 
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-10 10 20 30 

Figure 7: Numerical solution of the system ([2]) in case when the invariant kink-like 
solution dH]) with a = k = 1, mi = 1, m2 = 2, = 3, Ci = 1, C2 = 100, C3 = 0.01 
is taken as the Cauchy data. Successive graphs of TW solution, moving from left to 
right, correspond to tj = 3 (i — 1), i = I, ...5 




Figure 8: Numerical solution of the system ([2]) in case when the invariant kink-like 
solution ffTOj) with a = k = 1, mi = 0.25, m2 = 1713 = 1, Ci = C2 = C3 = 1 is taken 
as the Cauchy data. Successive graphs of TW solution, moving from left to right, 
correspond to = 5.5 {i — 1), i = 1, ...7 
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Figure 9: Numerical solution of the system ([2]) in case when the invariant kink-like 
solution (fTT]) with a = k = 1, m = 1, C2 = 1, C3 = 2 is taken as the Cauchy 
data. Successive graphs of TW solution, moving from left to right, correspond to 
ti = 3.75 (i- 1), i = 1,...7 

perturbation described by the formula f llOl) . Finally, the Fig. describes the numerical 
evolution of the N-shaped soliton with the heavy "tail", described by the formula (fTTl) . 
Results of numerical simulation show that the first two wave patterns evolve without 
any drastic changes of their shapes. In the last case diminishing of the maximal and 
minimal amplitudes is evidently seen. Besides the wave instability, this effect can be 
caused by the numerical scheme's viscosity. 

4 Summary 

So in this paper stability of TW solutions, satisfying the equation ([2]) under some 
restrictions on the values of the parameters are analyzed, and sufficient conditions for 
the stability of some family of exact solutions are presented in explicit form. It is 
quite evident, that the stability conditions stated by the theorem 2 are sufficient, but 
not necessary. In fact, they are the most strong among all possible conditions of this 
sort. No wonder, thus, that for some invariant TW solutions which do not satisfy 
the inequalities f l2^ . we succeeded to observe the stable self-similar evolution, as well. 
To gain the theoretical justification of the stability of the exact solutions differen from 
(128|) . an extra qualitative investigations, based on the more subtle methods are needed. 
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